!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==============================================================================================================================================|
recursive subroutine regrid (current)
!==============================================================================================================================================|

   use parameters
   use block_define
   use interfaces

   implicit none

   type(block), pointer, intent(inout) :: current

   integer(4) :: i, i_block, ip1, im1
   integer(4) :: flag

   real(8) :: dva(sys), dva2(sys), d2va2(sys), lv_exp, thral1, num_tmp
   type(block), pointer :: child

!----------------------------------------------------------------------------------------------------------------------------------------------|
   if (current%lv .lt. n_levels) then
      do i_block = 1, current%nc
         child => current%child(i_block)
         if (associated (child)) call regrid (child)
      enddo
   endif

   if (current%flag_skip .eq. 1) then
      current%pp%flag_skip = 1
      return
   endif

!----------------------------------------------------------------------------------------------------------------------------------------------|
!     error estimation (flag = 1 then block generation, flag = -1 then block destruction and flag = 0 for nothing)
!----------------------------------------------------------------------------------------------------------------------------------------------|
   flag = 0

!   if (current%lv .eq. 1) then
!      do i = ng + 1, current%nx - ng
!         ip1 = i + 1
!         im1 = i - 1
!         dro = abs (current%ro(ip1) - current%ro(im1)) / (abs (current%ro(i)) + eps)
!         dmx = abs (current%mx(ip1) - current%mx(im1)) / (abs (current%mx(i)) + eps)
!         dmy = abs (current%my(ip1) - current%my(im1)) / (abs (current%my(i)) + eps)
!         dmz = abs (current%mz(ip1) - current%mz(im1)) / (abs (current%mz(i)) + eps)
!         dbx = abs (current%bx(ip1) - current%bx(im1)) / (abs (current%bx(i)) + eps)
!         dby = abs (current%by(ip1) - current%by(im1)) / (abs (current%by(i)) + eps)
!         dbz = abs (current%bz(ip1) - current%bz(im1)) / (abs (current%bz(i)) + eps)
!         den = abs (current%en(ip1) - current%en(im1)) / (abs (current%en(i)) + eps)
!         if (dro .gt. thr_ro0 .or. dmx .gt. thr_mx0 .or. dmy .gt. thr_my0 .or. dmz .gt. thr_mz0 .or.                               &
!             dbx .gt. thr_bx0 .or. dby .gt. thr_by0 .or. dbz .gt. thr_bz0 .or. den .gt. thr_en0) then
!            flag = 1
!            len_regrid = min (len_regrid, min (current%x(i) - current%x(1), current%x(nx_block) - current%x(i)))
!            exit
!         else if (dro .gt. thr_ro0 * thr_de .or. dmx .gt. thr_mx0 * thr_de .or. dmy .gt. thr_my0 * thr_de .or.                     &
!             dmz .gt. thr_mz0 * thr_de .or. dbx .gt. thr_bx0 * thr_de .or. dby .gt. thr_by0 * thr_de .or.                          &
!             dbz .gt. thr_bz0 * thr_de .or. dbz .gt. thr_bz0 * thr_de) then
!            flag = 0
!         endif
!      enddo
!   else
!      parent => current%pp
!      do i = 1, nx_block
!         ro_p(i) = parent%ro(i)
!         mx_p(i) = parent%mx(i)
!         my_p(i) = parent%my(i)
!         mz_p(i) = parent%mz(i)
!         bx_p(i) = parent%bx(i)
!         by_p(i) = parent%by(i)
!         bz_p(i) = parent%bz(i)
!         en_p(i) = parent%en(i)
!      enddo
!      flux_bnd = parent%flux_bnd
!      select case (scheme)
!      case (1)
!         call tvd_mmc (ro_p, mx_p, my_p, mz_p, bx_p, by_p, bz_p, en_p, nx_block, ng,                                               &
!            flux_bnd, n_neighs, sys, parent%dx, parent%dt, gm, max_speed)
!      case (2)
!         call tvd_lf (ro_p, mx_p, my_p, mz_p, bx_p, by_p, bz_p, en_p, nx_block, ng,                                                &
!            flux_bnd, n_neighs, sys, parent%dx, parent%dt, gm, riemann_solver_flag)
!      case (3)
!         call weno (ro_p, mx_p, my_p, mz_p, bx_p, by_p, bz_p, en_p, nx_block, ng,                                                  &
!            flux_bnd, n_neighs, sys, parent%dx, parent%dt, gm, riemann_solver_flag)
!      endselect

!      do i = 1, nx_block
!         ro(i) = current%ro(i)
!         mx(i) = current%mx(i)
!         en(i) = current%en(i)
!      enddo
!      flux_bnd = current%flux_bnd
!      select case (scheme)
!      case (1)
!         call tvd_mmc (ro, mx, my, mz, bx, by, bz, en, nx_block, ng,                                                               &
!            flux_bnd, n_neighs, sys, current%dx, current%dt, gm, max_speed)
!      case (2)
!         call tvd_lf (ro, mx, my, mz, bx, by, bz, en, nx_block, ng,                                                                &
!            flux_bnd, n_neighs, sys, current%dx, current%dt, gm, riemann_solver_flag)
!      case (3)
!         call weno (ro, mx, my, mz, bx, by, bz, en, nx_block, ng,                                                                  &
!            flux_bnd, n_neighs, sys, current%dx, current%dt, gm, riemann_solver_flag)
!      endselect

!      lv_exp = 2 ** (thr_ex * parent%lv)
!      nx_h = nx_no_bnd / parent%nc
!      do i = 1, nx_h
!         i_p = i + nx_h * (current%p_pos - 1) + ng
!         i_c = i * 2 + ng
!         i_cm1 = i_c - 1
!         dro = abs (ro_p(i_p) - (ro(i_cm1) + ro(i_c)) / 2.0d0) / (abs (ro_p(i_p) + eps)) * lv_exp
!         dmx = abs (mx_p(i_p) - (mx(i_cm1) + mx(i_c)) / 2.0d0) / (abs (mx_p(i_p) + eps)) * lv_exp
!         dmy = abs (my_p(i_p) - (my(i_cm1) + my(i_c)) / 2.0d0) / (abs (my_p(i_p) + eps)) * lv_exp
!         dmz = abs (mz_p(i_p) - (mz(i_cm1) + mz(i_c)) / 2.0d0) / (abs (mz_p(i_p) + eps)) * lv_exp
!         dbx = abs (bx_p(i_p) - (bx(i_cm1) + bx(i_c)) / 2.0d0) / (abs (bx_p(i_p) + eps)) * lv_exp
!         dby = abs (by_p(i_p) - (by(i_cm1) + by(i_c)) / 2.0d0) / (abs (by_p(i_p) + eps)) * lv_exp
!         dbz = abs (bz_p(i_p) - (bz(i_cm1) + bz(i_c)) / 2.0d0) / (abs (bz_p(i_p) + eps)) * lv_exp
!         den = abs (en_p(i_p) - (en(i_cm1) + en(i_c)) / 2.0d0) / (abs (en_p(i_p) + eps)) * lv_exp
!         if (dro .gt. thr_ro .or. dmx .gt. thr_mx .or. dmy .gt. thr_my .or. dmz .gt. thr_mz .or.                                   &
!             dbx .gt. thr_bx .or. dby .gt. thr_by .or. dbz .gt. thr_bz .or. den .gt. thr_en) then
!            flag = 1
!            len_regrid = min (len_regrid, min (parent%x(i_p) - parent%x(1), parent%x(nx_block) - parent%x(i_p)))
!            exit
!         else if (dro .gt. thr_ro * thr_de .or. dmx .gt. thr_mx * thr_de .or. dmy .gt. thr_my * thr_de .or.                        &
!            dmz .gt. thr_mz * thr_de .or. dbx .gt. thr_bx * thr_de .or. dby .gt. thr_by * thr_de .or.                              &
!            dbz .gt. thr_bz * thr_de .or. den .gt. thr_en * thr_de) then
!            flag = 0
!         endif
!      enddo
!   endif

   if (current%lv .lt. n_levels) then

      flag = -1

      lv_exp = dble(ratio) ** (thr_ex * (current%lv - 1))
      num_tmp = 1.0d0 / 2.0d0
      thral1 = 1.0d0 - thr_al
      do i = 2, nx_block - 1
         ip1 = i + 1
         im1 = i - 1

         dva2(1) = num_tmp * abs(current%ro(ip1) - current%ro(im1))
         dva2(2) = num_tmp * abs(current%mx(ip1) - current%mx(im1))
         dva2(3) = num_tmp * abs(current%my(ip1) - current%my(im1))
         dva2(4) = num_tmp * abs(current%mz(ip1) - current%mz(im1))
         dva2(5) = num_tmp * abs(current%bx(ip1) - current%bx(im1))
         dva2(6) = num_tmp * abs(current%by(ip1) - current%by(im1))
         dva2(7) = num_tmp * abs(current%bz(ip1) - current%bz(im1))
         dva2(8) = num_tmp * abs(current%en(ip1) - current%en(im1))

         d2va2(1) = abs(current%ro(ip1) - 2.0d0 * current%ro(i) + current%ro(im1))
         d2va2(2) = abs(current%mx(ip1) - 2.0d0 * current%mx(i) + current%mx(im1))
         d2va2(3) = abs(current%my(ip1) - 2.0d0 * current%my(i) + current%my(im1))
         d2va2(4) = abs(current%mz(ip1) - 2.0d0 * current%mz(i) + current%mz(im1))
         d2va2(5) = abs(current%bx(ip1) - 2.0d0 * current%bx(i) + current%bx(im1))
         d2va2(6) = abs(current%by(ip1) - 2.0d0 * current%by(i) + current%by(im1))
         d2va2(7) = abs(current%bz(ip1) - 2.0d0 * current%bz(i) + current%bz(im1))
         d2va2(8) = abs(current%en(ip1) - 2.0d0 * current%en(i) + current%en(im1))

         dva(1) = (thr_al * dva2(1) / (abs(current%ro(i)) + eps) +                                                              &
                   thral1 * d2va2(1) / (dva2(1) + thr_fl * (abs(current%ro(i)) + eps))) * lv_exp
         dva(2) = (thr_al * dva2(2) / (abs(current%mx(i)) + eps) +                                                              &
                   thral1 * d2va2(2) / (dva2(2) + thr_fl * (abs(current%mx(i)) + eps))) * lv_exp
         dva(3) = (thr_al * dva2(3) / (abs(current%my(i)) + eps) +                                                              &
                   thral1 * d2va2(3) / (dva2(3) + thr_fl * (abs(current%my(i)) + eps))) * lv_exp
         dva(4) = (thr_al * dva2(4) / (abs(current%mz(i)) + eps) +                                                              &
                   thral1 * d2va2(4) / (dva2(4) + thr_fl * (abs(current%mz(i)) + eps))) * lv_exp
         dva(5) = (thr_al * dva2(5) / (abs(current%bx(i)) + eps) +                                                              &
                   thral1 * d2va2(5) / (dva2(5) + thr_fl * (abs(current%bx(i)) + eps))) * lv_exp
         dva(6) = (thr_al * dva2(6) / (abs(current%by(i)) + eps) +                                                              &
                   thral1 * d2va2(6) / (dva2(6) + thr_fl * (abs(current%by(i)) + eps))) * lv_exp
         dva(7) = (thr_al * dva2(7) / (abs(current%bz(i)) + eps) +                                                              &
                   thral1 * d2va2(7) / (dva2(7) + thr_fl * (abs(current%bz(i)) + eps))) * lv_exp
         dva(8) = (thr_al * dva2(8) / (abs(current%en(i)) + eps) +                                                              &
                   thral1 * d2va2(8) / (dva2(8) + thr_fl * (abs(current%en(i)) + eps))) * lv_exp

         if (any (dva .gt. thr_va0)) then
            flag = 1
            exit
         else if (any (dva .gt. thr_va0 * thr_de)) then
            flag = 0
         endif
      enddo
!      endif
   endif

!----------------------------------------------------------------------------------------------------------------------------------------------|
!     grid generation or destruction
!----------------------------------------------------------------------------------------------------------------------------------------------|
   if (flag .ne. -1) then
      current%flag_skip = 1
      current%pp%flag_skip = 1
   endif

   max_level = max(max_level, current%lv)
   if (flag .eq. 1) then
      if (.not. associated (current%child) .and. current%lv .lt. n_levels) call block_generate (current)
      if (current%lv .lt. n_levels) then
         do i_block = 1, current%nc
            child => current%child(i_block)
            call regrid (child)
         enddo
      endif
   elseif (flag .eq. -1) then
      if (current%neigh_lv(1) .le. current%lv + 1 .and. current%neigh_lv(2) .le. current%lv + 1 .and. associated (current%child)) then
         call block_destroy (current)
      endif
   endif

   return
end subroutine regrid
